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ABSTRACT 

This  paper  presents  a procedure  for  using  NASTRAN  to 
determine  the  flow  field  about  arbitrarily  shaped  bodies 
in  the  presence  of  a free  surface.  The  fundamental  un- 
known of  the  problem  is  the  velocity  potential  which  must 
satisfy  Laplace's  equation  in  the  fluid  region.  Boundary 
conditions  on  the  free  surface  may  involve  second  order 
derivatives  in  space  and  time.  In  cases  involving  infinite 
domains  either  a tractable  radiation  condition  is  applied 
at  a truncated  boundary  or  a series  expansion  is  used  and 
matched  to  the  local  finite  elements.  Solutions  are  pre- 
sented for  harmonic,  transient,  and  steady  state  problems 
and  compared  to  either  exact  solutions  or  other  numerical 
sol utions . 


INTRODUCTION 

The  pressure  distribution  and  flow  field  about  submerged  bodies  are 
important  in  the  determination  of  hydrodynamic  variables  such  as  lift  and  wave 
resistance  and  the  calculation  of  boundary-layer  characteristics.  The 
investigation  of  these  variables  can  be  realistically  modeled  by  assuming  the 
fluid  to  be  inviscid  and  incompressible.  In  this  case  the  equations  of  motion 
can  be  reduced  to  the  solution  of  Laplace's  equation  in  the  fluid  region.  The 
linearized  free  surface  condition  (small  wave  amplitude)  may  involve  second 
derivatives  of  the  velocity  potential  $ in  both  space  and  time  and  considerably 
complicates  the  problem.  The  free  surface  flows  investigated  in  this  paper  can 
be  divided  into  three  areas:  harmonic,  transient,  and  steady  state. 

An  exhaustive  list  of  literature  for  forced  harmonic  motion  or  diffraction 
problems  may  be  found  in  Wehausen  (ref.  1).  Problems  of  this  type  were 
generally  solved  by  using  a distribution  of  sources  or  dipoles  on  the  body 
boundary  with  an  appropriate  Green's  function  for  the  problem.  The  boundary 
condition  on  the  body  is  used  to  determine  the  strength  of  the  source  distri- 
bution (for  example,  Hess  and  Smith,  ref.  2).  Such  solutions  are  appropriate  only 
for  problems  of  infinite  or  constant  depth. 

Bai  (refs.  3-6)  uses  finite  elements  to  model  both  harmonic  and  steady  state 
problems  of  arbitrary  geometry.  Similar  methods  which  employ  variational 
functionals  have  been  used  by  Berkhoff  (ref.  7)  and  Chen  and  Mei  (ref.  8).  For 
steady  state  problems  Bai  developed  a localized  finite  element  method  (ref.  6) 
in  which  finite  elements  are  used  in  a localized  region  around  the  body  and  a 
series  expansion  is  used  in  the  remainder.  The  finite  element  representation 
is  matched  to  the  series  expansion  along  the  common  boundary  to  form  a consistent 
set  of  equations  for  the  nodal  potentials  and  series  coefficients. 

Finite  element  solutions  for  transient  free  surface  flows  are  given  by 
Visser  and  van  der  Wilt  (ref.  9).  Unfortunately,  for  the  transient  problem 
there  seems  to  be  no  suitable  method  to  construct  a completely  absorbing 
boundary  in  cases  involving  radiation  conditions.  For  that  reason  truncated 
boundaries  are  taken  far  enough  away  so  as  not  to  affect  the  region  of  interest. 
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The  purpose  of  this  paper  is  to  demonstrate  how  the  structural  analysis 
computer  program  NASTRAN  (refs.  10-11)  may  be  used  to  implement  finite  element 
procedures  for  modeling  the  three  types  of  free  surface  problems  described 
above.  The  use  of  NASTRAN  is  motivated  by  the  wide-ranging  capability,  con- 
venience of  use,  and  availability  of  this  general  purpose  computer  code.  The 
variety  of  finite  elements  available  in  NASTRAN  permits  the  method  presented 
here  for  2-D  and  axisymmetric  problems  to  be  routinely  applied  to  complex  3-D 
geometries  of  naval  and  marine  interest.  In  contrast  with  specialized  programming 
efforts,  NASTRAN  implementation  of  the  finite  element  procedures  is  enhanced  by 
a variety  of  pre-  and  postprocessing  programs  (ref.  12)  which  include  capabilities 
for  automatic  data  generation,  data  checking  via  interactive  graphics,  matrix 
bandwidth  and  profile  reduction  via  grid  point  resequencing,  and  contour  plots 
(in  the  case  of  scalar  variables  such  as  velocity  potential)  of  computer  output. 

In  addition,  the  NASTRAN  capability  to  model  free  surface  flow  problems  is 
currently  being  exploited  to  investigate  the  coupled  fluid-structure  interaction 
problem  involving  fluid  flow  about  an  elastic  body  near  or  on  a free  surface. 


FREE  SURFACE  EQUATIONS 


For  an  inviscid,  incompressible  fluid  in  an  irrotational  flow  field,  the 
equations  of  motion  and  continuity  reduce  to 

v2$  = 0 (1 ) 

where  ♦ is  the  velocity  potential  (ref.  13).  The  pressure  p in  the  fluid  can 
be  determined  from  Bernoulli's  equation. 


. E = 


p 


(2) 


where  p is  the  density  of  the  fluid  and  g is  the  gravitational  constant.  In 
Fig.  1 the  deflection  of  the  free  surface  n is  assumed  to  be  small  compared 
to  the  depth  d.  In  that  case  the  linearized  conditions  on  the  free  surface  are 
(ref.  10) 
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Once  the  potential  $ is  determined,  the  surface  elevation  n may  be  determined 
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from  Eq.  (4). 
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Figure  1.  Free  Surface  Wave 


HARMONIC  FREE  SURFACE  PROBLEMS 


2-D  Wave  Maker 

Consider  the  2-D  wavp  maker  shown  in  Fiq.  2.  At  x=0,  a wall  is  oscillating 
in  simple  harmonic  motion  with  velocity  V.  For  the  harmonic  problems,  assume 
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Figure  2.  Geometry  and  Boundary  Conditions  for  2-D  Wave  Maker 
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Then  Eqs.  (6)  and  (1)  qive 

V2$  = 0 

Usinq  Eqs.  (5),  (6),  and  (7),  the  free  surface  condition  becomes 


At  the  wall. 


and,  alonq  the  bottom. 
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The  solution  of  this  problem  can  be  obtained  by  separation  of  variables 
(see  Bai,  ref.  3)  and  is  qiven  by 


”a01't  ” ”aN* 

$(x.y)  = Aq  cosh  aQ(y+d)e  + AN  cos  aN(y+d)e 


‘O' 


where 


T 3 '*0  tanh(aNd) 


aN  tan(a^d)  for  all  N 


. - 4i  

0 sinh(2agd)  * 2.tgd 


'N  sin(2a^dT + 2a^d 


sinh(agd) 

a0 

sin(a^d) 

aN 


(12) 

(13) 

(14) 

(15) 

(16) 


The  first  term  of  Eq.  (12)  represents  a travelinq  wave  in  the  x-direction, 
while  the  succeedinq  terms  are  local  terms  that  are  only  siqnificant  for  small 
x.  Thus 
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for  large  x.  Eq.  (17)  is  a tractable  radiation  condition  which  can  be  applied 
at  suitable  boundary  far  enough  away  from  the  wall.  Eq.  (8),  together  with 
boundary  conditions  given  by  Eqs.  (9),  (10),  (11),  and  (17),  constitute  a well- 
posed  problem  for  Laplace's  equation. 

The  boundary  conditions  Eqs.  (9),  (10),  (11),  and  (17)  all  have  the  form 

+ Y*  = 6 (18) 

where  y and  s are  constants.  The  functional  form  for  Laplace's  equation  with 
the  mixed  boundary  condition  of  Eq.  (18)  is 

FU)  = hj\{^)2  + (ly)2| dA  + / {h*2~  B*)ds  (19) 

where  B is  the  boundary  of  region  A.  When  variations  are  taken  with  respect 
to  $ such  that 


6F  = 0 (20) 

then  Eqs.  (8)  and  (18)  are  satisfied. 


Eqs.  (19)  and  (20)  can  be  approximated  with  finite  elements  using  NASTRAN 
structural  elements.  A orocedure  for  using  structural  elements  to  model  fluid 
domains  which  satisfy  the  wave  equation  (or,  as  a special  case,  Laplace's 
equation)  is  given  by  Everstine  et  al  (ref.  14),  and  has  been  successfully 
applied  using  NASTRAN  on  several  problems  by  Schroeder  and  Marcus  (ref.  15), 
Marcus  (ref.  16),  and  Everstine  (ref.  17).  A translational  degree  of  freedom 
(in  this  case  the  x displacement)  is  chosen  to  represent  the  potential 
and  all  other  degrees  of  freedom  at  a node  are  permanently  constrained.  The 
linear  isoparametric  membrane  element,  QDMEM1  (NASTRAN  Level  16),  was  used. 

The  material  matrix  G and  the  mass  density  p of  the  QDMEM1  elements  are  chosen 
as  follows:  e 


G 


1 -1  0 
-1  1 0 

0 0 1 


p = 0 

e 


(21) 


NASTRAN' s Rigid  Format  8,  with  governing  equation  given  by 


(~u^M  + i o>B  + K)  £ = £(w) 


(22) 


is  chosen  as  the  analysis  method.  The  stiffness  matrix  K generated  by  the 
QDMEM1  elements  with  material  properties  given  by  Eq.  (21)  is  equivalent  to 
the  finite  element  representation  of  the  first  term  in  Eq.  (19). 
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The  free  surface  condition,  Eq.  (9), 
the  second  term  of  Eq.  (19).  A consistent 
implemented  using  NASTRAN  by  inserting  the 


(M2PP ) . . 

1 »J 


AX 

6g 


corresponds  to  y = (u2/g  and  B = 0 in 
formulation  for  this  term  is 
matrix 


1 

2 


i=k,t 

j=k,t 


(23) 


into  the  mass  matrix  M in  Eq.  (22)  using  DMIG  data  cards.  In  Eq.  (23),  k and 
i represent  the  two  nodes  which  lie  on  the  free  surface  for  each  of  the  QDMEM1 
surface  elements,  while  Ax  is  the  spacing  between  nodes  k and  l.  The  frequency 
u>  is  inserted  into  Eq.  (22)  using  a FREQ  data  card. 


The  radiation  condition,  Eq.  (17),  corresponds  to  b = 0 and  y = i ciq 
Eq.  (19).  A consistent  formulation  is  obtained  by  inserting  the  matrix  u 


(M2PP),  , = 


2 1 
1 2 


i =k , Jt 

j=M 


in 

(24) 


into  the  mass  matrix  M in  Eq.  (22)  using  DMIG  cards.  In  Eq.  (24)  k,  f,  and  Ay 
are  defined  as  in  Eq.  (23)  except  that  in  this  case  the  relevant  boundary 
surface  is  the  truncated  boundary. 


The  bottom  condition,  Eq.  (11),  is  a natural  boundary  condition  which  is 
automatically  satisfied  within  the  finite  element  approximation.  The  boundary 
condition  at  the  wall,  Eq.  (10),  is  implemented  by  inserting  the  vector 


Fi 


V Ay 


1/2 

1/2 


i=k,n 


(25) 


into  the  forcing  function  £(w)  in  Eq.  (22)  using  DAREA  data  cards.  The 
relevant  boundary  for  the  quantities  k,  i,  and  Ay  in  Eq.  (25)  is  the  oscillating 
wall . 


The  above  procedure  was  used  to  compute  the  fluid  response  for  the 
oscillating  wall  problem  illustrated  in  Fig.  2.  All  data  is  presented  in  non- 
dimensionalized  form  using  the  length  L and  the  velocity  V.  Results  are 
shown  in  Fig.  3 for  dimensionless  spacing  Ax  = Ay  = .0625  which  corresponds  to 
approximately  10  nodes  per  wave  length  for  the  linear  elements.  In  Fig.  3, 
the  amplitude  of  the  surface  elevation  is  plotted.  The  NASTRAN  solutions 
obtained  by  both  consistent  and  lumped  formulations,  as  well  as  the  analytic 
solution,  are  presented.  The  lumped  formulation  is  determined  by  using 
diagonalized  matrices  in  Eqs.  (23)  and  (24)  where  diagonal  terms  are  determined 
by  adding  together  all  terms  in  the  corresponding  row.  The  consistent  formula- 
tion is  a significant  improvement  over  the  lumped  formulation.  In  subsequent 
problems  only  a consistent  formulation  will  be  used. 
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Figure  3.  Amplitude  of  Surface  Elevation  for 
the  3-D  Wave  Maker 


Ax i symmetric  Wave  Maker 


The  geometry  and  boundary  conditions  for  the  axisymmetrit  wave  maker  ar 
shown  in  tig.  4.  Boundary  conditions  are  the  same  as  for  the  2-0  wave  maker 
except  for  the  additional  term  in  the  radiation  condition.  The  radiation 
condition  is  determined  by  investigating  the  exact  solution  (see  Bai,  ref.  3 


<f(r,z)  = B0Hn(,inr)  cosha0(zrd)  + >:  BN  HQ 


aN  i r)  cos  nN  (z+d)  ( 


4 sinh(a0d) 

B0  “ KjTa0r0)  sinh  [aQi)  + 3Sjj37 


where 


Fiqure  4.  Geometry  and  Boundary  Conditions  for 
Pulsating  Cylinder 


4i  sin(«Nd) 

BN  = ff[r-aNi  r0)aN{sin  aNd  + 2 aNd  (28) 

and  where  ag  and  are  given  by  Eqs.  (13)  and  (14),  rg  is  the  inner  radius  of 
the  cylinder,  and  Hg,  H-j  are  Hankel  functions  of  the  second  kind  of  order  0 
and  1,  respectively.  The  first  term  of  Eq.  (26)  is  an  outgoing  wave  and  the 
second  terms  represent  local  disturbances.  Thus  it  can  be  shown  that 

{ 2a  + i «0  1 4.  for  large  r (29) 

where  a is  defined  in  Fig.  4. 

This  problem  was  modeled  using  NASTRAN's  Rigid  Format  8.  CTRAPRG 
elements  were  used  (Everstine,  ref.  17)  with  dimensionless  spacing  given  by 
ax  = Ay  = .0625  (all  variables  are  non-dimensional ized  with  respect  to  V and  L). 
This  corresponds  to  approximately  10  nodes  per  wave  length.  Results  showing 
the  amplitude  of  the  surface  elevation  along  the  free  surface  are  presented  in 
Fig.  5.  These  results  are  based  on  applying  consistent  boundary  conditions, 
and  are  in  good  agreement  with  the  series  solution. 
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Figure  5.  Amplitude  of  Surface  Elevation  for 
the  Pulsating  Cylinder 


Refraction  Problems 


A surface  wave,  given  by 


n(x,t)  = A e 


i (wt-agX ) 


is  incident  upon  the  bottom  obstacle  shown  in  Fig.  6.  The  potential 
corresponding  to  the  incident  wave  is  given  by 


4’1(x,y)  = 


= A^l  C0Sh  a0  (y*d)  e 


-i  aQx 


COSh(ctgd) 


(30) 


(31) 


where  u>,  ag,  g.and  d satisfy  Eq.  (13).  In  order  to  determine  the  total 
potential  <j>  of  the  fluid  corresponding  to  the  incident  wave,  the  potential  <j> 
is  divided  into 

<t>  = + *R  (32) 
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Figure  6.  Geometry  and  Boundary  Conditions  for 
Refracted  Waves  Due  to  a Bottom  Obstacle 


where  4>R  is  the  refracted  potential.  The  boundary  conditions  and  governing 
eguations  on  are  shown  in  Fig.  6.  All  variables  are  non-dimens- onal i zed 
with  respect  to  the  length  L and  freguency  w,  and  boundary  conditions  are 
specified  in  a consistent  formulation. 

The  NASTRAN  results  shown  in  Figs.  7 and  8 are  presented  for  dimensionless 
spacing  Ax  = .125  and  Ay  = .0625, which  corresponds  to  approximately  41  nodes 
per  wave  length.  These  results  compare  favorably  with  the  finite  element 
solution  recently  re-computed  by  Bai  as  a correction  to  his  originally 
published  (ref.  5)  results.  Accuracies  within  4"  have  also  been  obtained  using 
coarser  grids  of  10-20  nodes  per  wave  length. 

A similar  free  surface  problem  is  illustrated  in  Fig.  9.  The  dimensionless 
spacing  used  was  Ax  = Ay  = .125, which  corresponds  to  approximately  42  nodes  per 
wave  length.  Again,  the  NASTRAN  results  shown  in  Figs.  10  and  11  compare  well 
with  the  finite  element  solution  recently  re-computed  by  Bai  (ref.  5). 


TRANSIENT  PROBLEMS 


Consider  the  transient  free  surface  problem  shown  in  Fig.  1?  illustrating 
the  time  dependent  pressure  distribution  on  the  free  surface,  The  pressure 
distribution  is  given  by 


p(x,t)  = P(x)  sin  ut 


(33) 
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Figure  10.  Amplitude  of  Refracted  Waves  for  the  Surface  Obstacle 

Eqs.  (35)  and  (36)  must  be  put  in  terms  of  $ and  3<t>/a t since  this  is  the  only 
suitable  input  to  NASTRAN.  Specifying  3n/3t  on  y=0  is  equivalent  by  Eq.  (3)  to 
specifying  3<t>/ ay  on  y=0.  Then  Laplace's  equation  may  be  solved  with  the 
boundary  conditions  shown  in  Fig.  12,  except  that  3$/3y  is  specified  on  the  free 
surface.  This  will  determine  » everywhere  initially.  Similarly,  specifying  n 
on  y=0  is  equivalent  by  Eq.  (4)  to  specifying  3$/3t  on  y=0.  Then  the  procedure 
just  described  may  be  repeated  to  determine  3$/3t  everywhere  initially,  since 
3$/3t  also  satisfies  Laplace's  equation  and  the  boundary  conditions  shown  in 
Fiq.  12  (not  including  the  free  surface  condition).  This  determines  ♦ and 
3<t>/3 1 everywhere  initially. 

The  variational  form  for  the  free  surface  problem  shown  in  Fig.  12,  based 
on  Hamilton’s  principle  (see  Courant  and  Hilbert,  ref.  18),  is 

^1  ^1  *1 

F(«)  ■ 1/  / I <5I>2«02M»df  / J(Lv»;*M>dsdf  / / i(|i)Jdxdt 

c 0 A l*X  OR  n <-9  3t 


+ / / ~r  TfT  ♦ dx  dt 
0 Free  pg  3t 
Surface 


0 Free  y 
Surface 
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Figure  11.  Phase  of  Refracted  Waves  for  the  Surface  Obstacle 

where  B is  the  boundary  of  the  region  A,  and  all  geometric  boundary  conditions 
are  enforced.  If  variations  of  F($)  taken  with  respect  to  * equal  zero. 

Eqs.  (1).  (5),  (18),  (35),  and  (36)  are  satisfied  for  zero  initial  conditions 
(f)*f2=0  in  Eqs.  (35)  and  (36)).  Non-zero  initial  conditions  can  be  easily 
incorporated  into  Eq.  (37). 

The  finite  element  representation  based  on  Eq.  (37)  was  implemented  using 
NASTRAN  by  modeling  the  fluid  with  QDMEM1  elements  where  material  properties 
are  given  by  Eq.  (21).  Any  translational  degree  of  freedom  can  be  used  to 
correspond  to  ♦,  but  all  remaining  degrees  of  freedom  are  permanently  con- 
strained. The  analysis  method  chosen  is  NASTRAVs  Rigid  Format  9,  with  the 
governing  equation  given  by 


Mi  + Bi  + kt  = F(t)  (38) 

The  stiffness  matrix  K generated  by  the  QDMEM1  elements  is  equivalent  to  the 
finite  element  representation  of  the  first  term  of  Eq.  (37). 

The  last  two  terms  of  Eq.  (37)  represent  the  free  surface  condition  and 
may  be  incorporated  into  NASTRAN  as  follows:  Let  t for  any  point  on  the  free 
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Figure  12.  Stationary  Pressure  Distribution  Oscillating  on  Free  Surface 
surface  be  yiven  by 


» ■ !'  N|  f.  (39) 

where  Ni  is  the  shape  function  for  node  i and  f.j  is  the  nodal  potential.  Then 
the  finite  element  formulation  for  the  third  term  of  Eg.  (37)  is  Implemented 
using  NASTRAN  by  inserting  the  matrix 

(M2PP)h  - ’ / N.N.  dx  (40) 

J 9 Free  J 
Surface 

into  the  mass  matrix  M in  Eg.  (38).  The  finite  element  representation  of  the 
last  term  of  Eg.  (37)  is  Implemented  using  NASTRAN  by  inserting  the  vector 

F.  • - / 1 $ N,dx  (41) 

] r 1 eg  >t  1 

Free 

Surface 

into  the  forcing  function  F(t)  In  Eg.  (38). 

Referring  to  Fig.  12,  the  natural  boundary  condition  a*/3n  » 0 
(corresponding  to  > 01  on  the  bottom  and  lett  face  is  automat icul 1> 
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Figure  13.  Development  of  Surface  Elevations  for  the 
Pressure  Distribution  Problem 


satisfied.  The  geometric  boundary  condition  * = 0 is  implemented  by  constraining 
=0  at  all  nodes  i on  the  downstream  boundary. 

The  above  procedure  was  used  to  model  the  geometry  and  boundary  conditions 
shown  in  Fig.  12.  All  variables  have  been  put  in  dimensionless  form  using  the 
pressure  Pg,  the  depth  L,  and  the  gravitational  constant  g. 
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This  procedure  was  used  to  model  the  geometry  and  boundary  conditions 
shown  in  Fig.  12.  QOMt M 1 elements  were  used  with  dimensionless  spacing 

* 0.1  to  0.5,  and  Ay-  0.25;  this  would  correspond  to  approximately  13  to  60 
nodes  per  wave  length,  where  the  wave  length  is  based  on  the  steady  state 
problem.  A dimensionless  time  step  of  At  3 . 1 was  used.  Rules  of  thumb  for 
estimating  spacing  and  time  steps  are  given  by  Visser  and  van  der  Wilt  (ref.  9). 
In  this  case  approximately  60  time  steps  per  period  of  the  forcing  function 
wore  used. 

In  Fig.  13.  the  NASTRAN  results  are  compared  to  a Fourier  series  solution 
given  by  Haussling  and  Van  Eseltine  (ref.  191.  The  wave  neights  are  in  good 
agreement  with  the  series  solution  and  illustrate  the  capability  of  NASTRAN  to 
model  transient  water  wave  problems. 


STIADY  SI  Alt  PROBLEMS 


Consider  the  steady  state  problem  shown  in  Fig.  14  where  a cylinder  of 
diameter  L is  moving  at  constant  velocity  U below  the  free  surface.  Steady  state 
solutions  are  sought  for  which  all  variables  are  independent  of  time  when 
referenced  to  a coordinate  system  moving  with  the  body,  that  is,  the  x-y 
coordinate  system  shown  in  Fig.  !4.  In  this  coordinate  system  it  can  be  shown 
that  the  potential  $ must  satisfy  Laplace's  equation,  and  the  free  surface 
condition  expressed  in  Eg.  (6)  becomes  (with  ps0  on  free  surface) 

B . Ul  l!±  m;‘1 


The  boundary  condition  on  the  rigid  cylinder  shown  in  Fig.  14  is 

» -U  coso  (43) 

where  t>  is  the  angle  between  the  x-direction  and  the  normal  to  the  body  directed 
out  of  the  fluid.  No  upstream  waves  are  allowed  and  the  Froude  number. 


is  such  that  downstream  waves  are  allowed  (see  Bai,  ref.  6).  Considerable 
effort  was  devoted  to  developing  tractable  radiation  conditions  for  the  up- 
stream and  downstream  boundaries,  resulting  in  the  conclusion  that  none  were 
possible.  For  this  reason  a series  expansion  is  used  in  the  regions  beyond  the 
upstream  and  downstream  truncated  boundaries  and  matched  (at  these  boundaries) 
to  the  finite  element  solution.  This  technique  was  developed  and  successfully 
applied  by  Bai  for  both  steady  state  problems  (ref.  6)  and  frequency  response 
problems  (ref.  4).  A similar  finite  element-series  expansion  technique  for  an 
acoustical  fluid  has  been  implemented  using  NASTRAN  by  Zarda  ^ret  .V' 


tan  a.d 
J 


U*  U‘ 

— a.  = tanh  .i-,d  , — a . - 

9 0 0 g j 

Upstream  from  the  body, 

N+l 

* = E B.f.  (46) 

j=l  J J 


The  sign  in  the  exponential  is  ( + ) for  the  upstream  boundary,  and  (-)  for  the 
downstream  boundary.  Furthermore,  N is  the  number  of  series  terms  chosen 
(the  same  number  is  assumed  upstream  and  downstream,  although  this  is  not 
necessary),  and  d is  the  depth.  Eqs.  (44)  and  (46)  satisfy  Laplace's  equation 
and  the  boundary  conditions  on  y=0  and  y=-d.  The  first  N terms  represent  local 
terms  that  decay  away  from  the  cylinder,  and  the  last  two  terms  in  Eq.  (44) 
represent  an  outgoing  downstream  wave;  no  such  waves  are  allowed  in  the  upstream 
expansion. 


Consider  the  variational  functional  given  by 


Free 

Surface 


/ U cos  i ds 
Body 


(47) 


- I * d> 

Upstream  x=-xL 


* / [^] 

Downstream  1 


“y  - fe] 


X — X r 


g "3xJTB 
x=xn 


y=0 


+ 


y=0 


where  points  A and  B and  boundaries  x^  and  xr  are  defined  in  Fig.  14,  and 
n is  the  normal  to  the  boundary  directed  out  of  the  fluid.  If  independent 
variations  of  F with  respect  to  ^a.  and  dre  set  edual  to  zero,  then 
Laplace's  equation  and  the  boundary  conditions  shown  in  Fig.  14  are  satisfied, 
and  3$/3n  is  continuous  on  the  upstream  and  downstream  boundaries.  No 
variations  of  the  bracketed  expressions  in  Eq.  (47)  are  allowed,  and  these 
expressions  can  be  evaluated  in  terms  of  the  series  coefficients  by  taking  the 
appropriate  derivatives  using  Eqs.  (44)  and  (46).  This  will  increase  the 
number  of  unknowns  by  the  number  (2N+4)  of  series  coefficients.  The  correspond- 
ing additional  equations  come  from  requiring  that  the  potential  $ is 
continuous  at  the  upstream  and  downstream  boundaries.  Let  the  finite  element 
representation  at  the  truncated  boundaries  be  given  by 

NN 

♦ = E N.  (48) 

i=l  1 1 


where  NN  is  the  number  of  nodes  on  the  truncated  boundary.  Then,  for  continuity 
of  $ on  the  downstream  boundary,  it  is  required  that 


Z M;  = Z A f 
i=l  k=1  K k 


on  x=xr 


and  on  the  upstream  boundary 


1 Mi  = E BiA 

i=l  k=l  K k 


on  x=  x, 


Eq.  (49)  is  multiplied  by  f.,  j=l  to  N+2,  and  integrated  from  -d  to  0.  This 
gives  a system  of  equations^ 


where 


i=£l  Vi  = k=!  Ak 


G.  . = / N.f . dy 
U >d  i J 


On  x=xr 


j=  1 to  N+2 


i = 1 to  NN 
j-=  1 to  N+3 


Hik  = / f A d>  J.k  = 1 to  N+3 


Eqs.  (51)  are  N+2  equations  involving  the  N+3  unknowns  Ak.  Multiplying  Eq.  (49) 
by  fjsj+3  and  integrating  from  -d  to  0 aoes  not  determine  an  independent  equation 
since  f^+2  is  proportional  to  f^+3  for  fixed  x. 

Multiplying  Eq.  (50)  by  f j , j=l  to  N+2,  and  integrating  from  -d  to  0 gives 


^ Vi  = ^ Vk 


x = x. 


j = 1 to  N+2 


Eqs.  (54)  are  N+2  equations  in  the  N+l  unknowns  Bk.  The  additional  equation, 
determined  by  multiplying  Eq.  (50)  by  ^+2'  corresponds  to  the  condition  that  no 
upstream  waves  are  allowed  (see  Bai , ref.  6).  Eqs.  (51)  and  (54)  give  the 
additional  2N+4  equations  involving  the  2N+4  unknowns  Aj  and  Bj. 

The  procedure  just  described  can  be  modeled  using  NASTRAN.  CIS2D8  elements 
are  used  to  model  the  fluid  (see  refs.  21  and  22).  These  second  order  iso- 
parametric elements  with  the  material  properties  given  by  Eq.  (21)  determine  a 
stiffness  matrix  equivalent  to  the  finite  element  representation  of  the  first 
term  of  Eq.  (47). 

The  second  term  of  Eq.  (47)  is  modeled  using  additional  CIS2D8  elements 
along  the  free  surface  as  shown  in  Fig.  14.  For  these  elements,  the  height  in 
the  y-direction  is  unity,  and  all  nodes  having  the  same  value  of  x are  con- 
strained to  move  together.  This  is  equivalent  to  having  1-D  isoparametric 
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elements  along  the  free  surface.  The  material  properties  for  these  elements 
are  given  by  Eq.  (21)  except  that  the  material  matrix  G is  multiplied  by  the 
constant  factor  (-U-’/g). 

The  third  term  of  Eq.  (47)  represents  a loading  term.  It  is  implemented 
using  NASTRAN  by  entering 


F = - / U cos  6 N . ds 

1 Body  1 


(55) 


as  nodal  forces,  where  N^  is  the  shape  function  for  node  i on  the  body. 


The  fourth  and  fifth  terms  of  Eq.  (47)  represent  coupling  terms  at  the 
upstream  and  downstream  boundaries.  Using  Eqs.  (44)  and  (46)  to  determine  the 
normal  derivatives,  the  finite  element  modeling  yields,  for  the  downstream 
boundary. 


(K2PP ) . . 


N+3 

E 

j = l 


0 


/ 

-d 


3f  • 

. _JL 

3X 


Ni  dy 


x=xR 


i=  1 to  NN 
j=  1 to  N+3 


(56) 


where  the  matrix  K2PP  is  added  to  the  stiffness  matrix.  In  order  to  implement 
this  condition,  N+3  scalar  unknowns  Aj  are  created  using  SPOINT  data  cards. 

Then  the  matrix  term  (K2PP)j  , in  Eq.  (56)  refers  to  node  i on  the  downstream 
boundary  and  to  the  SPOINT  representation  of  the  unknown  Aj.  Similarly,  for  the 
upstream  boundary 


(K2PP)  • 
* 


N+l  0 3f. 

E / — * 

j=l  -d 


3x 


Ni  dy 


i = 1 to  NN 
j = 1 to  N+l 


(57) 


For  the  last  two  terms  in  Eq.  (47),  the  finite  element  representation 
yields 


,2  3f; 


<K2PP)B,j  = ' 


X=Xr 


j = 1 to  N+3 


(58) 


m2 

<k2PP>a.j  ‘ T » 


x=x, 


j = 1 to  N+l 


(59) 


Eqs.  (51),  (54),  and  (56)  through  (59)  are  entered  into  NASTRAN  using  DM1G  cards 
and  complete  the  set  of  equations  to  solve  for  the  nodal  potentials  and  the 
upstream  and  downstream  series  coefficients.  NASTRAN's  Rigid  Format  1 (Static 
Analysis)  does  not  accept  DMIG  cards.  Therefore,  Rigid  Format  9 was  used  for 
one  time  step.  (Since  no  mass  or  damping  matrix  exists,  static  equilibrium  is 
reached  for  any  time  step.) 
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Computations,  with  all  quantities  being  non-dimensional ized  with  respect 
to  the  cylinder  diameter  L,  velocity  U,  and  fluid  density  p,  were  carried  out 
using  NASTRAN  for  the  grids  shown  in  Figs.  14  and  15.  Each  mesh  has  approxi- 
mately the  same  number  of  unknowns  since  the  series  solution  is  used  for 
| x 1 2 3.0  on  the  coarse  grid  and  for  | x | - 1.5  on  the  fine  grid.  Approxi- 
mately 9 and  17  nodes  per  wave  length  were  used  for  the  coarse  and  fine  grids, 
respectively. 

Wave  height  along  the  free  surface  is  plotted  in  Fig.  16.  Results  for  both 
the  coarse  and  fine  NASTRAN  grids  are  seen  to  compare  favorably  with  a solution 
obtained  by  Giesing  and  Smith  (ret.  23)  using  a distribution  of  sources.  The 
solutions  shown  here  all  satisfy  the  condition  that  no  upstream  waves  are 
allowed.  (In  this  case,  since  the  Froude  number  based  on  the  depth  is  less  than 
one,  downstream  waves  are  generated.) 

The  pressure  distribution  on  the  cylinder  may  be  determined  from  Bernoulli's 
equation.  Assuming  the  flow  about  the  cylinder  is  steady,  then,  in  the  x-y 
coordinate  system  that  is  moving  with  the  body,  Eq.  (2)  becomes 

p = - pu  i (^r'M^)?i  (6o) 

Fig.  17  illustrates  a plot  of  the  dimensionless  pressure  as  a function  of  the  x 
coordinate  on  the  surface  of  the  cylinder.  Results  are  shown  for  both  tne  fine 
and  coarse  grids  shown  in  Figs.  15  and  16.  The  discontinuities  of  the  curves 
occur  at  element  junctures  on  the  cylinder.  Although  the  potential  $ is 
necessarily  continuous,  3$/3x  and  S^/ay  are  not  necessarily  continuous  within 
the  finite  element  approximation , and  discontinuities  in  these  terms  are 
magnified  in  determining  the  pressure  in  Eq.  (60).  Also  shown  in  Fig.  17  is  a 
table  showing  computed  values  of  the  wave  resistance  and  lift  coefficients, 

Cp  and  Cp,  defined  by 


(pU2L)Cp  = 

- / Pdy 

(61) 

Body 

(pU-’L)C.  = 

/ Pdx 

(62) 

L 

Body 

The  values  of  Cp  and  C[  computed  using  NASTRAN  compare  favorably  with  those 
given  by  Giesing  and  Smith  (ref.  20). 


CONCLUSIONS 


The  problems  illustrated  here  demonstrate  the  capability  of  NASTRAN  to 
successfully  model  linearized  free  surface  flow  problems  for  harmonic,  transient, 
and  steady  state  cases.  Although  the  results  presented  here  are  for  arbitrary 
2-D  and  axisymmetric  geometries,  the  procedures  described  are  uirectly 
applicable  to  3-D  flow  problems  and  readily  extendable  to  the  coupled  problem  of 
fluid  flow  about  an  elastic  body. 


X COORDINATE  ON  SURFACE  OF  CYLINDER 


Figure  17.  Pressure  Distribution  on  the  Cylinder 


The  steady-state  flow  due  to  a cylinder  moving  below  the  free  surface 
was  computed  using  the  technique  of  coupling  finite  elements  with  a classical 
method  at  an  appropriate  common  boundary.  Finite  elements  are  used  to  model 
irregular  geometry  over  to  some  specified  regular  boundary,  and  classical 
solution  methods  are  used  beyond  this  boundary.  The  coupling  of  the  series 
solutions  to  the  finite  element  model  may  be  regarded  as  determining  a stiffness 
matrix  for  a "classical  finite  element."  Such  "elements",  if  available  in  the 
libraries  of  finite  element  computer  codes,  would  broaden  the  range  of  problems 
efficiently  handled  using  finite  elements.  Furthermore,  the  enhancement  of  the 
NASTRAN  capability  described  here  may  be  used  to  investigate  the  coupled 
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problem  of  fluid  flow  about  an  elastic  body  near  or  on  a free  surface.  In 
such  a case  both  the  structure  and  surrounding  fluid  would  be  modeled  using 
existing  NASTRAN  elements  and  would  be  coupled  at  the  fluid-structure 
interface. 
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